# Carica le librerie necessarie
library(here) # Gestione dei percorsi dei file
library(corrplot) # Visualizzazione matrice di correlazione
library(nortest) # Test di normalità
library(lmtest) # Test diagnostici per modelli lineari
library(car) # Test di collinearità e altri test per modelli lineari
library(plotly) # Visualizzazioni interattive
library(heatmaply) # Heatmap interattivo
library(ggheatmap) # Heatmap con ggplot2
library(gridExtra) # Organizzazione di grafici con grid
library(ggfortify) # Visualizzazioni grafiche per modelli statistici
library(olsrr) # Diagnostica e grafici per modelli di regressione
library(ggplot2) # Creazione di grafici con ggplot2
library(viridis) # Colormaps viridis
library(glmnet) # Modelli di regressione con penalizzazione lasso ed elastic net
library(highcharter) # Creazione di grafici interattivi con Highcharts
library(reshape2) # Manipolazione dati
library(RColorBrewer) # Colormaps
library(tidyr) # Manipolazione dati
library(dplyr) # Manipolazione dati
# Configura il percorso di root per la generazione di report
knitr::opts_knit$set(root.dir = here("0_Materiale"))
# Leggi il dataset da un file di testo
dataset <- read.delim("basketball_teams.txt")
# Definisci il primo e l'ultimo anno del range da considerare per lo studio
FIRST <- 1976
LAST <- 2011
# Filtra il dataset secondo le condizioni specificate
df <- dataset[dataset$lgID == "NBA" & dataset$year >= FIRST & dataset$year <= LAST & dataset$games == 82,]
# Converti la colonna lgID in un fattore per consentire la generazione di variabili dummy
dataset$lgID <- as.factor(dataset$lgID)
# Stampiamo un riassunto statistico del dataframe filtrato
summary(df)
## year lgID tmID franchID
## Min. :1976 Length:841 Length:841 Length:841
## 1st Qu.:1985 Class :character Class :character Class :character
## Median :1993 Mode :character Mode :character Mode :character
## Mean :1993
## 3rd Qu.:2001
## Max. :2008
## confID divID rank confRank
## Length:841 Length:841 Min. :0.000 Min. : 1.000
## Class :character Class :character 1st Qu.:2.000 1st Qu.: 4.000
## Mode :character Mode :character Median :3.000 Median : 7.000
## Mean :3.565 Mean : 7.164
## 3rd Qu.:5.000 3rd Qu.:10.000
## Max. :8.000 Max. :15.000
## playoff name o_fgm o_fga
## Length:841 Length:841 Min. :2565 Min. :5972
## Class :character Class :character 1st Qu.:2981 1st Qu.:6592
## Mode :character Mode :character Median :3220 Median :6903
## Mean :3239 Mean :6941
## 3rd Qu.:3489 3rd Qu.:7253
## Max. :3980 Max. :8868
## o_ftm o_fta o_3pm o_3pa o_oreb
## Min. :1189 Min. :1475 Min. : 0.0 Min. : 0.0 Min. : 720
## 1st Qu.:1523 1st Qu.:2039 1st Qu.: 78.0 1st Qu.: 293.0 1st Qu.: 974
## Median :1649 Median :2201 Median :283.0 Median : 814.0 Median :1083
## Mean :1666 Mean :2212 Mean :279.1 Mean : 804.6 Mean :1086
## 3rd Qu.:1797 3rd Qu.:2371 3rd Qu.:445.0 3rd Qu.:1267.0 3rd Qu.:1191
## Max. :2388 Max. :3051 Max. :837.0 Max. :2283.0 Max. :1520
## o_dreb o_reb o_asts o_pf o_stl
## Min. :2044 Min. :2922 Min. :1422 Min. :1476 Min. : 455.0
## 1st Qu.:2348 1st Qu.:3381 1st Qu.:1760 1st Qu.:1777 1st Qu.: 613.0
## Median :2433 Median :3506 Median :1934 Median :1900 Median : 671.0
## Mean :2434 Mean :3520 Mean :1934 Mean :1909 Mean : 681.3
## 3rd Qu.:2527 3rd Qu.:3647 3rd Qu.:2094 3rd Qu.:2033 3rd Qu.: 746.0
## Max. :2966 Max. :4216 Max. :2575 Max. :2470 Max. :1059.0
## o_to o_blk o_pts d_fgm d_fga
## Min. : 910 Min. :204.0 Min. : 6901 Min. :2488 Min. :5638
## 1st Qu.:1207 1st Qu.:360.0 1st Qu.: 7958 1st Qu.:2978 1st Qu.:6593
## Median :1311 Median :411.0 Median : 8404 Median :3243 Median :6911
## Mean :1338 Mean :421.8 Mean : 8423 Mean :3239 Mean :6941
## 3rd Qu.:1443 3rd Qu.:473.0 3rd Qu.: 8879 3rd Qu.:3493 3rd Qu.:7268
## Max. :2011 Max. :716.0 Max. :10371 Max. :4265 Max. :8142
## d_ftm d_fta d_3pm d_3pa d_oreb
## Min. :1217 Min. : 0.0 Min. : 0.0 Min. :1579 Min. : 745
## 1st Qu.:1514 1st Qu.: 82.0 1st Qu.: 282.0 1st Qu.:2033 1st Qu.: 986
## Median :1648 Median :280.0 Median : 836.0 Median :2203 Median :1095
## Mean :1666 Mean :279.1 Mean : 804.6 Mean :2212 Mean :1086
## 3rd Qu.:1808 3rd Qu.:446.0 3rd Qu.:1251.0 3rd Qu.:2395 3rd Qu.:1177
## Max. :2377 Max. :683.0 Max. :1768.0 Max. :3071 Max. :1495
## d_dreb d_reb d_asts d_pf d_stl
## Min. :2012 Min. :2976 Min. :1336 Min. :1434 Min. :461.0
## 1st Qu.:2326 1st Qu.:3378 1st Qu.:1778 1st Qu.:1788 1st Qu.:623.0
## Median :2431 Median :3497 Median :1939 Median :1900 Median :677.0
## Mean :2434 Mean :3520 Mean :1934 Mean :1909 Mean :681.3
## 3rd Qu.:2529 3rd Qu.:3653 3rd Qu.:2092 3rd Qu.:2020 3rd Qu.:734.0
## Max. :3067 Max. :4309 Max. :2537 Max. :2453 Max. :955.0
## d_to d_blk d_pts o_tmRebound d_tmRebound
## Min. : 949 Min. :264.0 Min. : 6909 Min. :0 Min. :0
## 1st Qu.:1208 1st Qu.:380.0 1st Qu.: 7968 1st Qu.:0 1st Qu.:0
## Median :1304 Median :419.0 Median : 8453 Median :0 Median :0
## Mean :1338 Mean :421.8 Mean : 8423 Mean :0 Mean :0
## 3rd Qu.:1444 3rd Qu.:460.0 3rd Qu.: 8841 3rd Qu.:0 3rd Qu.:0
## Max. :1980 Max. :654.0 Max. :10723 Max. :0 Max. :0
## homeWon homeLost awayWon awayLost neutWon
## Min. : 6.00 Min. : 1.00 Min. : 1.00 Min. : 8.00 Min. :0
## 1st Qu.:21.00 1st Qu.:10.00 1st Qu.:10.00 1st Qu.:21.00 1st Qu.:0
## Median :26.00 Median :15.00 Median :15.00 Median :26.00 Median :0
## Mean :25.63 Mean :15.37 Mean :15.37 Mean :25.63 Mean :0
## 3rd Qu.:31.00 3rd Qu.:20.00 3rd Qu.:20.00 3rd Qu.:31.00 3rd Qu.:0
## Max. :40.00 Max. :35.00 Max. :33.00 Max. :40.00 Max. :0
## neutLoss confWon confLoss divWon divLoss
## Min. :0 Min. : 5.00 Min. : 7.00 Min. : 1.00 Min. : 1.00
## 1st Qu.:0 1st Qu.:20.00 1st Qu.:20.00 1st Qu.: 8.00 1st Qu.: 9.00
## Median :0 Median :27.00 Median :26.00 Median :12.00 Median :12.00
## Mean :0 Mean :26.91 Mean :26.91 Mean :12.27 Mean :12.27
## 3rd Qu.:0 3rd Qu.:34.00 3rd Qu.:33.00 3rd Qu.:16.00 3rd Qu.:16.00
## Max. :0 Max. :48.00 Max. :52.00 Max. :25.00 Max. :27.00
## pace won lost games min
## Min. : 0.00 Min. :11 Min. :10 Min. :82 Min. :19680
## 1st Qu.: 0.00 1st Qu.:31 1st Qu.:32 1st Qu.:82 1st Qu.:19780
## Median : 0.00 Median :42 Median :40 Median :82 Median :19805
## Mean : 6.71 Mean :41 Mean :41 Mean :82 Mean :19817
## 3rd Qu.: 0.00 3rd Qu.:50 3rd Qu.:51 3rd Qu.:82 3rd Qu.:19855
## Max. :102.00 Max. :72 Max. :71 Max. :82 Max. :20080
## arena attendance bbtmID
## Length:841 Min. : 0 Length:841
## Class :character 1st Qu.:32767 Class :character
## Mode :character Median :32767 Mode :character
## Mean :32728
## 3rd Qu.:32767
## Max. :32767
# STUDIO DATI RELATIVI AI RIMBALZI
# o_oreb = Rimbalzi ottenuti in attacco
# o_dreb = Rimbalzi subiti in attacco
# o_reb = totale rimbalzi in attacco
# d_oreb = Rimbalzi subiti in difesa
# d_dreb = Rimbalzi ottenuti in difesa
# d_reb = totale rimbalzi in difesa
# Seleziona le colonne pertinenti dal dataframe
df_reb <- subset(df, select = c("o_oreb", "o_dreb", "o_reb", "d_oreb", "d_dreb", "d_reb", "won"))
# Imposta il layout a 2 righe e 3 colonne per i grafici di densità
par(mfrow = c(2, 3))
# Ciclo per generare grafici di densità per ciascuna variabile di interesse
for (variables in 1:(dim(df_reb)[2]-1)) {
thisvar <- df_reb[, variables]
d <- density(thisvar)
xmin <- floor(min(thisvar))
xmax <- ceiling(max(thisvar))
# Crea il plot della densità con stile accattivante
plot(d, main = names(df_reb)[variables], xlab = "", col = "blue", lwd = 1.5, xlim = c(xmin, xmax), ylim = c(0, max(d$y)*1.1))
# Aggiungi la distribuzione normale teorica ideale in rosso
x <- seq(xmin, xmax, length = 100)
lines(x, dnorm(x, mean = mean(thisvar), sd = sd(thisvar)), col = "red", lwd = 1.5) # Modifica lo spessore delle linee
# Aggiungi griglia per migliorare la leggibilità
grid()
}
# Aggiungi titolo in grassetto e corsivo, con spessore del testo modificato
title(bquote(bold(italic("Density plots with Normal Distribution"))), line = -17, cex.main = 2, outer = TRUE)
# Stampa percentuale di valori non zero per ciascuna variabile di interesse
print(paste("Percentage non-zero o_oreb: ", round(length(which(df_reb$o_oreb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero o_oreb: 100"
print(paste("Percentage non-zero o_dreb: ", round(length(which(df_reb$o_dreb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero o_dreb: 100"
print(paste("Percentage non-zero o_reb: ", round(length(which(df_reb$o_reb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero o_reb: 100"
print(paste("Percentage non-zero d_oreb: ", round(length(which(df_reb$d_oreb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero d_oreb: 100"
print(paste("Percentage non-zero d_dreb: ", round(length(which(df_reb$d_dreb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero d_dreb: 100"
print(paste("Percentage non-zero d_reb: ", round(length(which(df_reb$d_reb > 0)) / dim(df_reb)[1] * 100, 2)))
## [1] "Percentage non-zero d_reb: 100"
### SANDBOX
# Creazione di un vettore per archiviare un elenco di squadre
teams <- c()
# Creazione di una nuova variabile 'reb' come somma di 'o_reb' e 'd_reb'
df$reb <- c(df$o_reb + df$d_reb)
# Visualizzazione di un riepilogo statistico sulla variabile 'reb'
summary(df$reb)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5963 6805 7002 7040 7245 8359
# Loop per ottenere un elenco unico di squadre
for (team in df$tmID) {
teams <- unique(c(teams, team))
}
# Ordinamento delle squadre
teams <- sort(teams)
# Calcolo di valori aggregati per le variabili specificate
values <- aggregate(cbind(reb, o_oreb, o_dreb, d_oreb, d_dreb, o_reb, d_reb, won) ~ tmID, data = df, FUN = sum)
# Creazione di dati per una distribuzione normale
x <- seq(-4, 4, by = 0.01)
y <- dnorm(x, mean = 0, sd = 1)
normal <- data.frame(x, y)
# ISTOGRAMMA INTERATTIVO
histogram <- plot_ly(data = df, x = ~reb, type = "histogram",
marker = list(color = 'skyblue', line = list(color = 'black', width = 1)),
nbinsx = 20, legendgroup = "Rimbalzi", name = "Campione") %>%
add_trace(type = "scatter", mode = "lines",
x = c(mean(df$reb), mean(df$reb)),
y = c(0, 215),
line = list(color = "red", width = 2, dash = "dash"),
name = "Media") %>%
layout(title = list(text = "<b><i>Istogramma dei Rimbalzi</i></b>", pad = 10),
xaxis = list(title = '<b><i>Rimbalzi</i></b>'),
yaxis = list(title = '<b><i>Frequenza</i></b>'),
legend = list(title = "<b><i>Legenda</i></b>", tracegrouporder = "reversed"))
(histogram)
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
## Warning: 'scatter' objects don't have these attributes: 'nbinsx'
## Valid attributes include:
## 'cliponaxis', 'connectgaps', 'customdata', 'customdatasrc', 'dx', 'dy', 'error_x', 'error_y', 'fill', 'fillcolor', 'fillpattern', 'groupnorm', 'hoverinfo', 'hoverinfosrc', 'hoverlabel', 'hoveron', 'hovertemplate', 'hovertemplatesrc', 'hovertext', 'hovertextsrc', 'ids', 'idssrc', 'legendgroup', 'legendgrouptitle', 'legendrank', 'line', 'marker', 'meta', 'metasrc', 'mode', 'name', 'opacity', 'orientation', 'selected', 'selectedpoints', 'showlegend', 'stackgaps', 'stackgroup', 'stream', 'text', 'textfont', 'textposition', 'textpositionsrc', 'textsrc', 'texttemplate', 'texttemplatesrc', 'transforms', 'type', 'uid', 'uirevision', 'unselected', 'visible', 'x', 'x0', 'xaxis', 'xcalendar', 'xhoverformat', 'xperiod', 'xperiod0', 'xperiodalignment', 'xsrc', 'y', 'y0', 'yaxis', 'ycalendar', 'yhoverformat', 'yperiod', 'yperiod0', 'yperiodalignment', 'ysrc', 'key', 'set', 'frame', 'transforms', '_isNestedKey', '_isSimpleKey', '_isGraticule', '_bbox'
# PLOT DI DENSITA CON SOVRAPPOSIZIONE DI UNA NORMALE IDEALE
density_data <- density(df$reb)
mu <- mean(df$reb)
sigma <- sd(df$reb)
normal_data <- dnorm(density_data$x, mean = mu, sd = sigma)
density_plot <- plot_ly(x = density_data$x, y = density_data$y, type = 'scatter', mode = 'lines',
line = list(color = 'blue', width = 2),
name = "Densità rimbalzi totale") %>%
layout(title = "Densità dei Rimbalzi",
xaxis = list(title = "Rimbalzi"),
yaxis = list(title = "Density", autotick = TRUE, autorange = TRUE))
density_plot <- add_trace(density_plot, x = density_data$x, y = normal_data, mode = 'lines',
line = list(color = 'green', width = 2),
fill = "tozeroy", fillcolor = "rgba(0, 255, 0, 0.2)",
name = "Dist normale ideale")
density_plot <- add_trace(density_plot, x = c(mu, mu), y = c(0, max(density_data$y)),
mode = 'lines', line = list(color = 'red', width = 2, dash = 'dash'),
name = "Media")
(density_plot)
# CORRPLOT
data_subset <- df[, c("reb", "o_reb", "d_reb")]
cor_matrix <- cor(data_subset)
rownames(cor_matrix) <- colnames(data_subset)
colnames(cor_matrix) <- colnames(data_subset)
cor_data <- reshape2::melt(cor_matrix)
names(cor_data) <- c("Var1", "Var2", "Corr")
dimnames(cor_matrix) <- list(rownames(cor_matrix), colnames(cor_matrix))
data_for_plotly <- as.data.frame(as.table(cor_matrix))
cor_plot <- plot_ly(data = cor_data,
x = ~Var1,
y = ~Var2,
z = ~Corr,
type = "heatmap",
colors = colorRampPalette(c("#4575b4", "#91bfdb", "#e0f3f8", "#fee08b", "#d73027"))(100),
hoverinfo = "x+y+z") %>%
layout(title = 'Correlation Matrix',
xaxis = list(title = "", tickangle = 45, side = "bottom", automargin = TRUE),
yaxis = list(title = "", automargin = TRUE),
autosize = TRUE)
cor_values <- round(as.matrix(cor_matrix), 2) # Round for readability
for (i in seq_len(nrow(cor_matrix))) {
for (j in seq_len(ncol(cor_matrix))) {
cor_plot <- cor_plot %>% add_annotations(
x = rownames(cor_matrix)[i],
y = colnames(cor_matrix)[j],
text = as.character(cor_values[i, j]),
showarrow = FALSE,
font = list(color = ifelse(cor_values[i, j] < 0.5, "white", "black"))
)
}
}
(cor_plot)
# BOX PLOT
aggregated_data <- df %>%
select(tmID, year, o_reb, d_reb) %>%
group_by(tmID, year) %>%
summarize(
o_reb = mean(o_reb, na.rm = TRUE),
d_reb = mean(d_reb, na.rm = TRUE),
.groups = "drop" # Aggiunto per evitare il raggruppamento
)
reshaped_data <- aggregated_data %>%
pivot_longer(cols = c(o_reb, d_reb), names_to = "stat_type", values_to = "stat") %>%
unite("new_col", tmID, stat_type, sep = "_") %>%
pivot_wider(names_from = new_col, values_from = "stat")
box_plot <- plot_ly()
for (i in 2:ncol(reshaped_data)) {
current_column_data <- reshaped_data[[i]]
box_plot <- box_plot %>% add_trace(y = current_column_data, name = colnames(reshaped_data)[i], type = "box")
}
(box_plot)
## Warning: Ignoring 30 observations
## Warning: Ignoring 30 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 4 observations
## Warning: Ignoring 4 observations
## Warning: Ignoring 23 observations
## Warning: Ignoring 23 observations
## Warning: Ignoring 8 observations
## Warning: Ignoring 8 observations
## Warning: Ignoring 24 observations
## Warning: Ignoring 24 observations
## Warning: Ignoring 12 observations
## Warning: Ignoring 12 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 29 observations
## Warning: Ignoring 29 observations
## Warning: Ignoring 30 observations
## Warning: Ignoring 30 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 9 observations
## Warning: Ignoring 9 observations
## Warning: Ignoring 26 observations
## Warning: Ignoring 26 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 3 observations
## Warning: Ignoring 3 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 21 observations
## Warning: Ignoring 21 observations
## Warning: Ignoring 11 observations
## Warning: Ignoring 11 observations
# HEATMAP
heatmap_df <- subset(values, select = -c(o_oreb, o_dreb, d_oreb, d_dreb, won, reb))
rownames(heatmap_df) <- heatmap_df$tmID
heatmap_df <- heatmap_df[,-1]
heatmap_plot <- heatmaply(
heatmap_df,
colors = viridis(n = 256, option = "magma"),
k_col = 2,
k_row = 4,
)
(heatmap_plot)
# BARPLOT
bar_plot <- plot_ly(values, x = ~tmID, y = ~o_reb, type = 'bar', name = 'Rimbalzi offensivi', marker = list(color = '#FFAFA1')) %>%
add_trace(y = ~d_reb, name = 'Rimbalzi difensivi', marker = list(color = '#b2fff8')) %>%
layout(yaxis = list(title = 'Valori'), barmode = 'stack')
(bar_plot)
#TEST ANDERSON-DARLING
# Questo codice esegue il test di Anderson-Darling sulla variabile 'reb' nel tuo dataset (df)
ad.test(df$reb)
##
## Anderson-Darling normality test
##
## data: df$reb
## A = 3.6997, p-value = 3.1e-09
#Con un livello di significatività (α) di 0.01 e un p-value molto piccolo (3.1e-09) ottenuto dal test di normalità di Anderson-Darling per i dati della variabile df\(reb, puoi concludere che hai sufficiente evidenza statistica per respingere lipotesi nulla che i dati seguono una distribuzione normale.Con il tuo livello di significatività del 0.01 e il p-value molto piccolo (3.1e-09), il p-value è inferiore al livello di significatività, quindi respingeresti lipotesi nulla. Questo suggerisce che i dati nella variabile df\)reb non seguono una distribuzione normale al livello di significatività del 0.01. In termini più pratici, hai abbastanza evidenza statistica per concludere che la variabile df$reb non segue una distribuzione normale basandoti sui risultati del test di Anderson-Darling.
#TEST KOLMOGOROV SMIRNOV
# Test di Kolmogorov-Smirnov per confrontare la distribuzione di 'reb' con una distribuzione normale
ks.test(df$reb, "pnorm")
## Warning in ks.test.default(df$reb, "pnorm"): i legami non dovrebbero essere
## presenti per il test di Kolmogorov-Smirnov
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: df$reb
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided
#TEST SHAPIRO WILK
# Test di Shapiro-Wilk per la normalità dei dati nella variabile 'reb' (rimbalzi)
sf.test(df$reb)
##
## Shapiro-Francia normality test
##
## data: df$reb
## W = 0.98016, p-value = 1.758e-08
#In sintesi, il risultato del test di Shapiro-Francia indica che i tuoi dati nella variabile df$reb non seguono una distribuzione normale. Questo è supportato dal valore basso del p-value, il quale suggerisce che la differenza tra la distribuzione dei tuoi dati e una distribuzione normale è statisticamente significativa.
#STUDIO SULLE VITTORIE
# ISTOGRAMMA INTERATTIVO
# Calcola la media della variabile "won" nel data frame "df"
mean_value <- mean(df$won)
# Crea un istogramma con plot_ly
histogram <- plot_ly(df, x = ~won, type = "histogram",
marker = list(color = "skyblue", line = list(color = "white", width = 0.5)),
opacity = 0.7, name = "Campione") %>%
layout(title = list(text = "<b><i>Distribuzione delle Vittorie</i></b>", y = 0.97),
xaxis = list(title = "<b><i>Numero di Vittorie</i></b>", zeroline = FALSE),
yaxis = list(title = "<b><i>Frequenza</i></b>", zeroline = FALSE),
barmode = "overlay") %>%
add_trace(x = ~mean(won), type = "scatter", mode = "lines",
line = list(color = "red", width = 2), name = "<i><b>Media</i></b>") %>%
add_trace(x = ~mean(won), type = "scatter", mode = "markers",
marker = list(color = "red", size = 8), showlegend = FALSE) %>%
add_annotations(text = sprintf("<i><b>Media: %.2f</b></i>", mean(df$won)), x = mean_value, y = 0,
arrowhead = 2, arrowcolor = "red", arrowsize = 1.5, arrowwidth = 2)
# Visualizza l'istogramma interattivo
histogram
## A marker object has been specified, but markers is not in the mode
## Adding markers to the mode...
# PLOT DENSITA INTERATTIVO
# Calcola la densità delle vittorie
density_plot <- density(df$won)
# Crea il plot di densità con plot_ly
density_interactive <- plot_ly(x = density_plot$x, y = density_plot$y, type = "scatter", mode = "lines",
line = list(color = "skyblue", width = 2), name = "Densità") %>%
layout(title = list(text = "<b><i>Distribuzione di Densità delle Vittorie</i></b>", x = 0.5),
xaxis = list(title = "<i>Numero di Vittorie</i>"),
yaxis = list(title = "<i>Densità</i>"),
showlegend = TRUE) %>%
add_trace(x = c(mean(df$won), mean(df$won)), y = c(0, max(density_plot$y)),
type = "scatter", mode = "lines", line = list(color = "red", width = 2),
name = "<i><b>Media</b></i>") %>%
add_trace(x = mean(df$won), y = max(density_plot$y), type = "scatter", mode = "markers",
marker = list(color = "red", size = 8), showlegend = FALSE) %>%
add_annotations(text = sprintf("<i><b>Media: %.2f</b></i>", mean(df$won)), x = mean(df$won), y = max(density_plot$y) * 1.05,
arrowhead = 2, arrowcolor = "red", arrowsize = 1.5, arrowwidth = 2,
ax = 0, ay = -40)
# Visualizza il plot di densità interattivo
density_interactive
## A line object has been specified, but lines is not in the mode
## Adding lines to the mode...
# CORRPLOT
data_subset <- df[, c(11:40, 54)]
cor_matrix <- cor(data_subset, use = "complete.obs")
cor_data <- melt(cor_matrix)
names(cor_data) <- c("Variable1", "Variable2", "Correlation")
# Crea la heatmap
p <- plot_ly(data = cor_data, x = ~Variable1, y = ~Variable2, z = ~Correlation,
type = "heatmap",
colors = viridis(n = 1024, option = "magma"),
hoverinfo = "x+y+z") %>%
layout(
title = '<b><i>Correlation Matrix Heatmap</i></b>',
xaxis = list(
title = list(text = "<b><i>Variabile 1</i></b>", standoff = 0),
tickangle = 45,
zeroline = FALSE
),
yaxis = list(
title = list(text = "<b><i>Variabile 2</i></b>", standoff = 0),
tickangle = 45,
zeroline = FALSE
),
autosize = TRUE
)
# Visualizza la heatmap
(p)
# BOXPLOT
aggregated_data <- df %>%
select(tmID, year, won) %>%
group_by(tmID, year) %>%
summarize(
won = mean(won, na.rm = TRUE),
)
## `summarise()` has grouped output by 'tmID'. You can override using the
## `.groups` argument.
reshaped_data <- aggregated_data %>%
pivot_longer(cols = c(won), names_to = "stat_type", values_to = "stat") %>%
unite("new_col", tmID, stat_type, sep = "_") %>%
pivot_wider(names_from = new_col, values_from = "stat")
fig <- plot_ly()
for (i in 2:ncol(reshaped_data)) {
current_column_data <- reshaped_data[[i]]
fig <- fig %>% add_trace(y = current_column_data, name = colnames(reshaped_data)[i], type = "box")
}
# Visualizza il boxplot interattivo
(fig)
## Warning: Ignoring 30 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 4 observations
## Warning: Ignoring 23 observations
## Warning: Ignoring 8 observations
## Warning: Ignoring 24 observations
## Warning: Ignoring 12 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 29 observations
## Warning: Ignoring 30 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 31 observations
## Warning: Ignoring 13 observations
## Warning: Ignoring 9 observations
## Warning: Ignoring 26 observations
## Warning: Ignoring 1 observations
## Warning: Ignoring 19 observations
## Warning: Ignoring 3 observations
## Warning: Ignoring 27 observations
## Warning: Ignoring 21 observations
## Warning: Ignoring 11 observations
df$reb <- df$o_reb + df$d_reb
# Questo codice effettua alcune operazioni di preprocessing sui dati, come la creazione di nuove variabili (f1 a f10) e la divisione del dataset in set di training e test. Successivamente, viene creato un modello lineare (linMod) e ne viene eseguito il resampling. Infine, viene creato un modello lineare con le covariate normalizzate (linModNormalized) e vengono generati grafici di diagnostica per entrambi i modelli.
# o_oreb = Rimbalzi ottenuti in attacco
# o_dreb = Rimbalzi subiti in attacco
# o_reb = totale rimbalzi in attacco = o_oreb + o_dreb
# d_oreb = Rimbalzi subiti in difesa
# d_dreb = Rimbalzi ottenuti in difesa
# d_reb = totale rimbalzi in difesa
# Definizione di nuove variabili
df$f1 <- (df$o_oreb) / (df$o_fga - df$o_fgm)
df$f2 <- (df$d_dreb) / (df$d_fga - df$d_fgm)
df$f3 <- (df$o_oreb + 1.5 * df$d_dreb) / (df$o_dreb + 2 * df$d_oreb)
df$f4 <- (df$o_oreb - df$o_dreb) + 1.5 * (df$d_dreb - df$d_oreb)
df$f5 <- (df$d_oreb / df$d_to) / (df$o_dreb / df$o_to)
df$f6 <- (df$o_oreb + df$d_dreb) - (df$d_oreb - df$o_dreb)^2
df$f7 <- (df$f1 + df$f2)^2
df$f8 <- (df$o_oreb) / (df$o_reb)
df$f9 <- ((df$o_oreb) / (df$o_dreb))^2
df$f10 <- ((df$d_dreb) / (df$d_oreb))^2
# Divisione in Test e Train per evitare che il modello fitti troppo bene sui nostri dati
sample <- sample(c(TRUE, FALSE), nrow(df), replace=TRUE, prob=c(0.7, 0.3))
df = subset(df, select = c("f1", "f2", "f3", "f4", "f5", "f6", "f7", "f8", "f9", "f10", "won", "divID", "confID"))
train <- df[sample, ]
test <- df[!sample, ]
#ATTENZIONE: facendo il train sul valore delle variabili, questo significa che sono esse ad essere i nostri dati, non i rimbalzi in se.
# importante non confondersi durante l'esposizione.
# Creazione del modello lineare
linMod <- lm(won ~ f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 + f10, data = train)
summary(linMod)
##
## Call:
## lm(formula = won ~ f1 + f2 + f3 + f4 + f5 + f6 + f7 + f8 + f9 +
## f10, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7696 -3.1466 0.0134 3.1314 13.3878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.021e+03 8.602e+01 11.874 < 2e-16 ***
## f1 -6.650e+02 1.481e+02 -4.491 8.52e-06 ***
## f2 -1.055e+03 1.565e+02 -6.744 3.65e-11 ***
## f3 -3.961e+02 3.630e+01 -10.911 < 2e-16 ***
## f4 7.341e-02 7.724e-03 9.504 < 2e-16 ***
## f5 -1.819e+02 4.754e+00 -38.253 < 2e-16 ***
## f6 1.236e-05 1.484e-06 8.328 5.64e-16 ***
## f7 4.622e+02 7.823e+01 5.908 5.81e-09 ***
## f8 -2.537e+02 7.397e+01 -3.429 0.000647 ***
## f9 3.281e+01 2.236e+01 1.467 0.142912
## f10 4.732e+00 9.798e-01 4.829 1.74e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.683 on 597 degrees of freedom
## Multiple R-squared: 0.8602, Adjusted R-squared: 0.8578
## F-statistic: 367.2 on 10 and 597 DF, p-value: < 2.2e-16
# Normalizziamo le covariate
train$f1_z <- scale(train$f1)
train$f2_z <- scale(train$f2)
train$f3_z <- scale(train$f3)
train$f4_z <- scale(train$f4)
train$f5_z <- scale(train$f5)
train$f6_z <- scale(train$f6)
train$f7_z <- scale(train$f7)
train$f8_z <- scale(train$f8)
train$f9_z <- scale(train$f9)
train$f10_z <- scale(train$f10)
# Normalizzo i dati di test
test$f1_z <- scale(test$f1)
test$f2_z <- scale(test$f2)
test$f3_z <- scale(test$f3)
test$f4_z <- scale(test$f4)
test$f5_z <- scale(test$f5)
test$f6_z <- scale(test$f6)
test$f7_z <- scale(test$f7)
test$f8_z <- scale(test$f8)
test$f9_z <- scale(test$f9)
test$f10_z <- scale(test$f10)
linModNormalized <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z + f7_z + f8_z + f9_z + f10_z, data = train)
# Imposta il layout della pagina
par(mfrow = c(2, 2))
# Grafico dei residui standardizzati
plot(linModNormalized, which = 1, col = "skyblue", pch = 16, main = "Residui Standardizzati")
abline(h = 0, col = "red", lty = 2) # Aggiungi una linea orizzontale a zero nel grafico dei residui standardizzati
# Grafico dei livelli
plot(linModNormalized, which = 2, col = "lightgreen", pch = 16, main = "Grafico dei Livelli")
# Grafico della distribuzione dei residui
plot(linModNormalized, which = 3, col = "pink", pch = 16, main = "Distribuzione dei Residui")
# Grafico di Q-Q plot dei residui
plot(linModNormalized, which = 4, col = "orange", pch = 16, main = "Q-Q Plot dei Residui")
# Aggiungi altre personalizzazioni a piacere, come titoli, etichette degli assi, ecc.
# Ripristina il layout predefinito
par(mfrow = c(1, 1))
# TEST SUL MODELLO DI REGRESSIONE LINEARE
# 1 Summary
summary(linModNormalized)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z +
## f7_z + f8_z + f9_z + f10_z, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7696 -3.1466 0.0134 3.1314 13.3878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.4030 0.1899 212.755 < 2e-16 ***
## f1_z -21.0213 4.6809 -4.491 8.52e-06 ***
## f2_z -50.2841 7.4564 -6.744 3.65e-11 ***
## f3_z -25.4430 2.3318 -10.911 < 2e-16 ***
## f4_z 19.5615 2.0582 9.504 < 2e-16 ***
## f5_z -12.8531 0.3360 -38.253 < 2e-16 ***
## f6_z 6.4973 0.7802 8.328 5.64e-16 ***
## f7_z 38.0074 6.4329 5.908 5.81e-09 ***
## f8_z -8.2180 2.3964 -3.429 0.000647 ***
## f9_z 2.0413 1.3915 1.467 0.142912
## f10_z 6.8481 1.4180 4.829 1.74e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.683 on 597 degrees of freedom
## Multiple R-squared: 0.8602, Adjusted R-squared: 0.8578
## F-statistic: 367.2 on 10 and 597 DF, p-value: < 2.2e-16
# 2 R-quadrato e R-quadrato Adattato
summary_linModNormalized <- summary(linModNormalized)
r_squared <- summary_linModNormalized$r.squared
cat("R-squared:", r_squared, "\n")
## R-squared: 0.8601693
n <- length(df$o_oreb)
k <- length(linModNormalized$coefficients) - 1
adjusted_r_squared <- 1 - ((1 - r_squared) * (n - 1) / (n - k - 1))
cat("Adjusted R-squared:", adjusted_r_squared, "\n")
## Adjusted R-squared: 0.9872881
# 2 Test Shapiro per valutare la normalità dei residui
shapiro.test(residuals(linModNormalized))
##
## Shapiro-Wilk normality test
##
## data: residuals(linModNormalized)
## W = 0.99432, p-value = 0.02275
# 3 Test di omoschedasticità (Breusch-Pagan test) --> il risultato suggerisce omoschedasiticità
bptest(linModNormalized)
##
## studentized Breusch-Pagan test
##
## data: linModNormalized
## BP = 14.173, df = 10, p-value = 0.1652
# 4 Test di multicollinearità
car::vif(linModNormalized)
## f1_z f2_z f3_z f4_z f5_z f6_z
## 606.573045 1539.138319 150.527628 117.266989 3.125347 16.850589
## f7_z f8_z f9_z f10_z
## 1145.588930 158.982557 53.601265 55.664014
#TROVO GLI OUTLAYERS
# Il codice implementa l'individuazione degli outliers considerando i residui che superano una soglia moltiplicata per la deviazione standard.
# Calcola i residui dal modello di regressione lineare normalizzato
residui <- residuals(linModNormalized)
# Definisci una soglia per gli outlier
soglia_outlier <- 2
outliers <- which(abs(residui) > soglia_outlier*sd(residui))
# Identifica gli outlier basati sulla soglia definita
outliers <- which(abs(residui) > soglia_outlier * sd(residui))
outliers
## 440 470 519 652 769 799 822 837 891 900 914 929 936 940 945 957
## 1 23 62 162 239 260 277 290 328 334 344 357 361 365 369 379
## 1026 1107 1128 1139 1205 1214 1224 1301
## 405 462 480 489 533 538 545 603
#MODELLO SENZA OUTLAYERS
#Il codice rimuove gli outliers, ricrea il modello lineare normalizzato e produce grafici di diagnostica per entrambi i modelli.
# Rimuovi gli outliers e ricrea il modello lineare normalizzato
if (length(outliers) != 0) {
train_1 <- train[-outliers,]
} else {
train_1 <- train
}
# ATTENZIONE: facendo il train sul valore delle variabili, questo significa che sono esse ad essere i nostri dati, non i rimbalzi in sé. Importante non confondersi durante l'esposizione.
linModNormalized_1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z + f7_z + f8_z + f9_z + f10_z, data = train_1)
# Stampa i summary dei modelli
summary(linModNormalized)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z +
## f7_z + f8_z + f9_z + f10_z, data = train)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7696 -3.1466 0.0134 3.1314 13.3878
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.4030 0.1899 212.755 < 2e-16 ***
## f1_z -21.0213 4.6809 -4.491 8.52e-06 ***
## f2_z -50.2841 7.4564 -6.744 3.65e-11 ***
## f3_z -25.4430 2.3318 -10.911 < 2e-16 ***
## f4_z 19.5615 2.0582 9.504 < 2e-16 ***
## f5_z -12.8531 0.3360 -38.253 < 2e-16 ***
## f6_z 6.4973 0.7802 8.328 5.64e-16 ***
## f7_z 38.0074 6.4329 5.908 5.81e-09 ***
## f8_z -8.2180 2.3964 -3.429 0.000647 ***
## f9_z 2.0413 1.3915 1.467 0.142912
## f10_z 6.8481 1.4180 4.829 1.74e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.683 on 597 degrees of freedom
## Multiple R-squared: 0.8602, Adjusted R-squared: 0.8578
## F-statistic: 367.2 on 10 and 597 DF, p-value: < 2.2e-16
summary(linModNormalized_1)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z +
## f7_z + f8_z + f9_z + f10_z, data = train_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0156 -2.9560 0.1958 3.0363 9.3808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.1581 0.1742 230.573 < 2e-16 ***
## f1_z -18.9076 4.7109 -4.014 6.77e-05 ***
## f2_z -48.9506 7.4952 -6.531 1.44e-10 ***
## f3_z -27.7666 2.1909 -12.674 < 2e-16 ***
## f4_z 22.5459 1.9693 11.449 < 2e-16 ***
## f5_z -13.0032 0.3099 -41.956 < 2e-16 ***
## f6_z 7.1655 0.7371 9.721 < 2e-16 ***
## f7_z 36.0699 6.4990 5.550 4.37e-08 ***
## f8_z -9.8409 2.1973 -4.479 9.07e-06 ***
## f9_z 2.5699 1.2670 2.028 0.043 *
## f10_z 7.4600 1.3101 5.694 1.98e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.206 on 573 degrees of freedom
## Multiple R-squared: 0.884, Adjusted R-squared: 0.882
## F-statistic: 436.7 on 10 and 573 DF, p-value: < 2.2e-16
# Grafici
ols_plot_resid_lev(linModNormalized)
ols_plot_resid_lev(linModNormalized_1)
# Influence plot per il primo modello
influencePlot(linModNormalized, id = 5)
## Warning in applyDefaults(id, defaults = list(method = "noteworthy", n = 2, :
## unnamed id arguments, will be ignored
## StudRes Hat CookD
## 769 -2.337140 0.140832493 0.080791911
## 822 2.567135 0.141797087 0.098069596
## 940 2.887305 0.007421224 0.005597557
## 957 3.037847 0.108873451 0.101105957
# Personalizzazione del primo grafico
par(mar=c(5, 5, 2, 2))
title(main = "Influence Plot - Modello 1", col.main = "blue", font.main = 4)
# Apri una nuova finestra grafica
dev.new()
# Influence plot per il secondo modello
influencePlot(linModNormalized_1, id = 5)
## Warning in applyDefaults(id, defaults = list(method = "noteworthy", n = 2, :
## unnamed id arguments, will be ignored
## StudRes Hat CookD
## 454 -1.92755178 0.07523360 2.734929e-02
## 592 0.41907733 0.12016398 2.183703e-03
## 643 2.26269646 0.04316141 2.084517e-02
## 903 -0.02277664 0.13549723 7.404726e-06
## 907 2.25139453 0.01156714 5.354468e-03
## 1206 1.56841793 0.10113810 2.509849e-02
# Personalizzazione del secondo grafico
par(mar=c(5, 5, 2, 2))
# Aggiungi un titolo sopra al grafico
title(main = "Influence Plot - Modello 2", col.main = "blue", font.main = 4)
# METODO DI REGRESSIONE LASSO
# Il codice implementa il metodo di regressione LASSO con cross-validation per trovare il miglior valore lambda.
# Definizione della variabile di risposta
y <- train_1$won
# Definizione della matrice delle variabili predittive (uso solo poche variabili ma potete farlo con tutte da togliere però la risposta area)
x <- data.matrix(train_1[, c("f1_z", "f2_z", "f3_z", "f4_z", "f5_z", "f6_z", "f7_z", "f8_z", "f9_z", "f10_z")])
# Esegui la cross-validation k-fold per trovare il valore lambda ottimale
cv_model <- cv.glmnet(x, y, alpha = 1)
# Trova il valore lambda ottimale che minimizza il test MSE
best_lambda <- cv_model$lambda.min
best_lambda
## [1] 0.003616601
# Addestramento del modello con cross-validation
cv_model <- cv.glmnet(x, y)
# Grafico personalizzato
plot(cv_model, xvar="lambda", main="", xlab="log(Lambda)", ylab="Mean Squared Error", col="blue", lwd=2)
## Warning in plot.window(...): parametro grafico "xvar" non valido
## Warning in plot.xy(xy, type, ...): parametro grafico "xvar" non valido
## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xvar" non valido
## Warning in axis(side = side, at = at, labels = labels, ...): parametro grafico
## "xvar" non valido
## Warning in box(...): parametro grafico "xvar" non valido
## Warning in title(...): parametro grafico "xvar" non valido
# Aggiungi un titolo personalizzato con formattazione LaTeX per grassetto e corsivo
title(main=expression(bold(italic("Validazione Incrociata"))), line = 3)
# Fittiamo il modello con il miglior lambda (penalizzazione)
best_model <- glmnet(x, y, alpha = 1, lambda = best_lambda)
# Stampa i coefficienti del modello LASSO
coef(best_model)
## 11 x 1 sparse Matrix of class "dgCMatrix"
## s0
## (Intercept) 40.139921
## f1_z -3.850296
## f2_z -24.088259
## f3_z -24.780290
## f4_z 19.280082
## f5_z -13.008981
## f6_z 6.166835
## f7_z 14.584547
## f8_z -10.168556
## f9_z 3.876559
## f10_z 5.714035
# Stampa il summary del modello lineare normalizzato precedente per confronto
summary(linModNormalized_1)
##
## Call:
## lm(formula = won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z +
## f7_z + f8_z + f9_z + f10_z, data = train_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.0156 -2.9560 0.1958 3.0363 9.3808
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 40.1581 0.1742 230.573 < 2e-16 ***
## f1_z -18.9076 4.7109 -4.014 6.77e-05 ***
## f2_z -48.9506 7.4952 -6.531 1.44e-10 ***
## f3_z -27.7666 2.1909 -12.674 < 2e-16 ***
## f4_z 22.5459 1.9693 11.449 < 2e-16 ***
## f5_z -13.0032 0.3099 -41.956 < 2e-16 ***
## f6_z 7.1655 0.7371 9.721 < 2e-16 ***
## f7_z 36.0699 6.4990 5.550 4.37e-08 ***
## f8_z -9.8409 2.1973 -4.479 9.07e-06 ***
## f9_z 2.5699 1.2670 2.028 0.043 *
## f10_z 7.4600 1.3101 5.694 1.98e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.206 on 573 degrees of freedom
## Multiple R-squared: 0.884, Adjusted R-squared: 0.882
## F-statistic: 436.7 on 10 and 573 DF, p-value: < 2.2e-16
# Questo codice esegue stime usando un modello Lasso e un modello lineare (LM) su un set di dati di test e ne valuta le prestazioni utilizzando l'errore quadratico medio (RMSE) e le metriche di classificazione binaria, ad esempio matrice di confusione, accuratezza, precisione, sensibilità, punteggio F e specificità.
# Predizioni utilizzando il modello Lasso
new <- data.matrix(test[, c("f1_z", "f2_z", "f3_z", "f4_z", "f5_z", "f6_z", "f7_z", "f8_z", "f9_z", "f10_z")])
prevLasso <- predict(best_model, s = best_lambda, newx = new)
# Predizioni utilizzando il modello LM (Linear Model)
new <- subset(test, select = c("f1_z", "f2_z", "f3_z", "f4_z", "f5_z", "f6_z", "f7_z", "f8_z", "f9_z", "f10_z"))
prevLM <- predict(linModNormalized_1, newdata = new)
# Calcolo dell'Errore Quadratico Medio (RMSE) per il modello Lasso
rmsLasso <- sqrt(mean((test$won - prevLasso)^2))
# Calcolo dell'Errore Quadratico Medio (RMSE) per il modello LM (Linear Model)
rmsLM <- sqrt(mean((test$won - prevLM)^2))
# Classificazione binaria utilizzando una soglia (0.5) per le predizioni del modello Lasso
prev <- ifelse(prevLasso > 0.5, "1", "0")
prev <- as.factor(prev)
# Matrice di Confusione per le predizioni del modello Lasso
confMatrix <- table(prev, test$won)
# Metriche di Performance
accuracy <- sum(confMatrix[1], confMatrix[4]) / sum(confMatrix[1:4])
precision <- confMatrix[4] / sum(confMatrix[4], confMatrix[2])
sensitivity <- confMatrix[4] / sum(confMatrix[4], confMatrix[3])
fscore <- (2 * (sensitivity * precision)) / (sensitivity + precision)
specificity <- confMatrix[1] / sum(confMatrix[1], confMatrix[2])
# Visualizza i Risultati
print(paste("Accuratezza:", round(accuracy, digits = 4)))
## [1] "Accuratezza: 0.6"
print(paste("Precisione:", round(precision, digits = 4)))
## [1] "Precisione: 0.6667"
print(paste("Sensibilità:", round(sensitivity, digits = 4)))
## [1] "Sensibilità: 0.6667"
print(paste("F-Score:", round(fscore, digits = 4)))
## [1] "F-Score: 0.6667"
print(paste("Specificità:", round(specificity, digits = 4)))
## [1] "Specificità: 0.5"
# Confronto tra RMS
print(paste("RMS Lasso:", round(rmsLasso, digits = 4)))
## [1] "RMS Lasso: 5.1466"
print(paste("RMS Modello Lineare:", round(rmsLM, digits = 4)))
## [1] "RMS Modello Lineare: 5.2047"
#TEST ANOVA
# Esecuzione di un test ANOVA per la variabile 'confID' come predittore
# Supponiamo che 'won' sia la variabile dipendente
resp_conf <- anova(lm(won ~ confID, data = train_1))
# Verifica se la variabile 'confID' ha un effetto significativo sulle vittorie
if (resp_conf["confID", "Pr(>F)"] < 0.05) {
print("Si rifiuta l'ipotesi nulla perché la variabile Conference ha un effetto sulle vittorie. Posso inserire a modello la variabile.")
} else {
print("Si accetta l'ipotesi nulla perché la variabile Conference non ha un effetto sulle vittorie. Non posso inserire a modello la variabile.")
}
## [1] "Si accetta l'ipotesi nulla perché la variabile Conference non ha un effetto sulle vittorie. Non posso inserire a modello la variabile."
# Esecuzione di un test ANOVA per la variabile 'divID' come predittore
# Supponiamo che 'won' sia la variabile dipendente
resp_div <- anova(lm(won ~ divID, data = train_1))
# Verifica se la variabile 'divID' ha un effetto significativo sulle vittorie
if (resp_div["divID", "Pr(>F)"] < 0.05) {
print("Si rifiuta l'ipotesi nulla perché la variabile Division ha un effetto sulle vittorie. Posso inserire a modello la variabile.")
div <- TRUE
} else {
print("Si accetta l'ipotesi nulla perché la variabile Division non ha un effetto sulle vittorie. Non posso inserire a modello la variabile.")
div <- FALSE
}
## [1] "Si accetta l'ipotesi nulla perché la variabile Division non ha un effetto sulle vittorie. Non posso inserire a modello la variabile."
#EFFETTI INTERAZIONE
# Verifica del valore di 'div' per costruire il modello appropriato
if (div) {
linModNormalized1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z + f7_z + f8_z + f9_z + f10_z + divID + f1_z:f2_z + f4_z:divID, data = train_1)
} else {
linModNormalized1 <- lm(won ~ f1_z + f2_z + f3_z + f4_z + f5_z + f6_z + f7_z + f8_z + f9_z + f10_z + f1_z:f2_z + f4_z:divID, data = train_1)
}
# Estrazione del riassunto del modello
riassunto <- summary(linModNormalized1)
# Estrazione dei nomi delle righe dal riassunto del modello
nomi_righe <- rownames(riassunto$coefficients)
# Nomi delle righe da escludere
nomi_da_escludere <- c("(Intercept)", "f1_z:f2_z", "f4_z:divIDCD", "f4_z:divIDMW", "f4_z:divIDNW", "f4_z:divIDPC", "f4_z:divIDSE", "f4_z:divIDSW", "divIDCD", "divIDMW", "divIDNW", "divIDPC", "divIDSE", "divIDSW")
# Creazione di un vettore con i nomi delle righe da utilizzare
nomi_righe_da_utilizzare <- setdiff(nomi_righe, nomi_da_escludere)
# Aggiunta della variabile 'divID' se presente nei nomi delle righe
if ("divIDCD" %in% nomi_righe) {
nomi_righe_da_utilizzare <- append(nomi_righe_da_utilizzare, c("divID"))
}
# Indici di inizio e fine per le p-values della variabile 'divID'
div_indice_inizio <- which(rownames(riassunto$coefficients) == "f4_z:divIDCD")
div_indice_fine <- which(rownames(riassunto$coefficients) == "f4_z:divIDSW")
# Estrazione delle p-values della variabile 'divID'
div_p_values <- riassunto$coefficients[div_indice_inizio:div_indice_fine, "Pr(>|t|)"]
# Aggiunta di 'f4_z:divID' se la media delle p-values è inferiore a 0.05
if (mean(div_p_values) < 0.05) {
variabili <- append(nomi_righe_da_utilizzare, c("f4_z:divID"))
}
# Aggiunta di 'f1_z:f2_z' se la sua p-value è inferiore a 0.05
if (riassunto$coefficients["f1_z:f2_z", "Pr(>|t|)"] < 0.05) {
variabili <- append(nomi_righe_da_utilizzare, c("f1_z:f2_z"))
}
# Se nessuna delle condizioni precedenti è soddisfatta, utilizza solo le variabili esistenti
if (!(mean(div_p_values) < 0.05) && !(riassunto$coefficients["f1_z:f2_z", "Pr(>|t|)"] < 0.05)) {
variabili <- nomi_righe_da_utilizzare
}
# Creazione della formula significativa per il modello
formula_significativa <- as.formula(paste("won ~", paste(variabili, collapse = " + ")))
# Costruzione del nuovo modello lineare con la formula significativa
linModNormalized1 <- lm(formula_significativa, data = train_1)
# Visualizzazione del riassunto del nuovo modello
summary(linModNormalized1)
##
## Call:
## lm(formula = formula_significativa, data = train_1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.6428 -2.8848 0.1359 2.9845 9.4655
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.9445 0.1972 202.572 < 2e-16 ***
## f1_z -21.3191 4.8113 -4.431 1.12e-05 ***
## f2_z -51.5660 7.5554 -6.825 2.25e-11 ***
## f3_z -26.5996 2.2420 -11.864 < 2e-16 ***
## f4_z 21.6454 2.0014 10.815 < 2e-16 ***
## f5_z -13.0055 0.3088 -42.117 < 2e-16 ***
## f6_z 6.6184 0.7726 8.566 < 2e-16 ***
## f7_z 38.5394 6.5652 5.870 7.38e-09 ***
## f8_z -9.3337 2.2005 -4.242 2.59e-05 ***
## f9_z 2.5617 1.2624 2.029 0.0429 *
## f10_z 6.4646 1.3764 4.697 3.31e-06 ***
## f1_z:f2_z -0.4396 0.1927 -2.281 0.0229 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.191 on 572 degrees of freedom
## Multiple R-squared: 0.8851, Adjusted R-squared: 0.8829
## F-statistic: 400.4 on 11 and 572 DF, p-value: < 2.2e-16
# Questo codice addestra un modello di regressione generalizzata di Poisson (linModNormalized2_pois) utilizzando la stessa specifica del modello lineare (linModNormalized2). Successivamente, vengono ottenuti e combinati i coefficienti di entrambi i modelli in un unico dataframe per una facile comparazione.
# Creazioni di un modello di regressione generalizzata di Poisson
linModNormalized1_pois = glm(formula_significativa, family=poisson(link=log), data = train_1)
summary(linModNormalized1_pois)
##
## Call:
## glm(formula = formula_significativa, family = poisson(link = log),
## data = train_1)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.6487423 0.0076641 476.082 < 2e-16 ***
## f1_z -0.0759700 0.1871492 -0.406 0.684792
## f2_z -0.6559339 0.2949496 -2.224 0.026156 *
## f3_z -0.8218722 0.0881466 -9.324 < 2e-16 ***
## f4_z 0.6768258 0.0787680 8.593 < 2e-16 ***
## f5_z -0.3390626 0.0120410 -28.159 < 2e-16 ***
## f6_z 0.1840628 0.0281955 6.528 6.66e-11 ***
## f7_z 0.3814903 0.2581319 1.478 0.139437
## f8_z -0.2773063 0.0821204 -3.377 0.000733 ***
## f9_z 0.0838793 0.0479279 1.750 0.080098 .
## f10_z 0.2164135 0.0520268 4.160 3.19e-05 ***
## f1_z:f2_z -0.0005862 0.0072869 -0.080 0.935880
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2294.39 on 583 degrees of freedom
## Residual deviance: 326.91 on 572 degrees of freedom
## AIC: 3552
##
## Number of Fisher Scoring iterations: 4
# Ottenimento dei coefficienti per ambo i modelli
normal = coefficients(linModNormalized1)
poisson = exp(coefficients(linModNormalized1_pois))
# Combinazione dei coefficienti in un unico dataframe
coefficients_table <- cbind(normal, poisson)
coefficients_table
## normal poisson
## (Intercept) 39.9445260 38.4263085
## f1_z -21.3191022 0.9268440
## f2_z -51.5659536 0.5189572
## f3_z -26.5996116 0.4396079
## f4_z 21.6454468 1.9676222
## f5_z -13.0055452 0.7124379
## f6_z 6.6184349 1.2020914
## f7_z 38.5393686 1.4644655
## f8_z -9.3336530 0.7578224
## f9_z 2.5617476 1.0874976
## f10_z 6.4645993 1.2416157
## f1_z:f2_z -0.4395592 0.9994139
#Considerazioni Generali:
# È importante notare che gli effetti delle variabili possono essere interpretati in modo diverso a seconda della distribuzione scelta per il modello. La scelta tra distribuzione normale e di Poisson dipende dalla natura della tua variabile dipendente e dai tuoi obiettivi di modellazione. Nel contesto di modelli di regressione, è sempre buona pratica verificare l'adeguatezza del modello esaminando i residui, eseguendo test diagnostici e valutando la bontà di adattamento. L'interpretazione dei coefficienti dovrebbe essere fatta considerando la scala appropriata per la distribuzione utilizzata (lineare per la normale, logaritmica per la Poisson). Se stai cercando di prevedere il numero di vittorie, la distribuzione di Poisson potrebbe essere più
# appropriata per variabili conteggio come questa. Tuttavia, è sempre necessario verificare l'adeguatezza del modello ai dati specifici.